function [output,solution] = produce_moments(Xi,Rho,Nu,to_plot,fignum)

%save('runs/local_save.mat');

% Xi            are estimated structural parameters
% Rho           are fixed structural parameters
% Nu            are discretization, estimation, and simulation parameters

% output        is the simulated output
% solution      is the model solution corresponding to [Xi,Rho,Nu]
   
% Reset random number generator
rng(Nu.globalrandseed);

%%%%%%%%%%%%%
% Parameters
%%%%%%%%%%%%%

% Assign estimated parameters
for p=1:Nu.Npars
        pars.(Nu.targeted_pars{p})                      = Xi(p);
end

% Assign fixed parameters
for p=1:numel(Nu.fixed_pars)
        pars.(Nu.fixed_pars{p})                         = Rho(p);
end

% If xi, the reduced-form parameter combining k and sigma, is estimated,
% update the value of k correspondingly
if (~any(strcmp(Nu.targeted_pars,'k')))&&( any(strcmp(Nu.targeted_pars,'xi')) || any(strcmp(Nu.targeted_pars,'sigma')) )
        pars.k                                          = (pars.sigma/pars.xi-1)^(-1)*pars.theta - pars.r;
end
if (any(strcmp(Nu.targeted_pars,'k')))
        pars.xi                                         = (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma;
end
if isequal(Nu.targeted_pars,{'S','sigma','M_e','C'})
        pars.xi                                         = (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma;
end
if ~( abs( (1+pars.theta/(pars.r+pars.k))^(-1)*pars.sigma - pars.xi ) < eps )
    save('temp/produce_moments_call.mat')
    error('Inconsistent parameters in produce_moments')
end

%%%%%%%%%%%%%%%%%
% Model solution
%%%%%%%%%%%%%%%%%

loud                                        = 0;
solution                                    = solve(pars,Nu,loud);

%%%%%%%%%%%%%%%%%%%
% Model simulation
%%%%%%%%%%%%%%%%%%%

Nt                                          = ceil(pars.tSMM/Nu.Deltat);
idx_shock                                   = ceil((5+Nu.shock_date/30)/Nu.Deltat);
tvec                                        = 0:Nu.Deltat:((Nt)*Nu.Deltat);
[tvec_eom,it_eom]                           = deal(NaN(pars.tSMM,1));
tvec_eom(1)                                 = 0;
it_eom(1)                                   = 1;
for t=1:(pars.tSMM-1)
        [~,it_eom(t+1)]                             = find(tvec<=t,1,'last');
        tvec_eom(t+1)                               = tvec(it_eom(t+1));
end
[simM,simX]                                 = deal(NaN(Nt+1,Nu.Nsim,Nu.Nreps));
[M_monthly,X_monthly]                       = deal(NaN(numel(it_eom),Nu.Nsim,Nu.Nreps));

% Initial conditions for X
InitX                                       = Nu.InitX;

% Moments
Moms                                        = NaN(numel(Nu.log_targeted_moments),Nu.Nreps);
%var_Moms                                    = NaN(8,8,Nu.Nreps);

% Simulate M_t
if Nu.mc_method == 0

    a                                           = [1 -(1-pars.theta*Nu.Deltat)];
    b                                           = pars.sigma*sqrt(Nu.Deltat);
    Phi                                         = solution.preT.Phi_0;
    Xvec                                        = solution.preT.Nu.X';
    
    % Simulate from the model
    for nr=1:Nu.Nreps

            epsilonSim                              = Nu.eps{nr};
            if isfield(Nu,'which_moments') && strcmp(Nu.which_moments,'baseline')
            epsilonSim(idx_shock,:)                 = (-pars.M_c*pars.S/4 + (pars.sigma_d*Nu.exposure)')/(pars.sigma*sqrt(Nu.Deltat));
            elseif isfield(Nu,'which_moments') && strcmp(Nu.which_moments,'persist')
            epsilonSim(idx_shock,:)                 = (-pars.M_c*pars.S/4 + (pars.sigma_d*Nu.exposure)')/(pars.sigma*sqrt(Nu.Deltat));
            else
            epsilonSim(idx_shock,:)                 = (-pars.M_c*pars.S + (pars.sigma_d*Nu.exposure)')/(pars.sigma*sqrt(Nu.Deltat));
            end
            localsimM                               = pars.M_c + filter(b,a,epsilonSim,0,1);
            simM(:,:,nr)                            = localsimM; 
            dummy                                   = 1;
            localstruct                             = simulate_X( pars, Nu, InitX, Nt, localsimM, Phi, Xvec, dummy);
            localsimX                               = localstruct.X;
            simX(:,:,nr)                            = localsimX;

           
            % Sample monthly
            M_monthly(:,:,nr)                       = localsimM(it_eom,:);
            localsimX_monthly                       = localsimX(it_eom,:);
            X_monthly(:,:,nr)                       = localsimX_monthly;

            % Moments 
            if isfield(Nu,'which_moments') && strcmp(Nu.which_moments,'baseline')
                Moms(:,nr)                              = compute_data_moments_mod(localsimX_monthly,Nu);
            elseif isfield(Nu,'which_moments') && strcmp(Nu.which_moments,'persist')
                Moms(:,nr)                              = compute_data_moments_IRF_mod_cf_persist(localsimX_monthly,Nu);
            else
                Moms(:,nr)                              = compute_data_moments_IRF_mod(localsimX_monthly,Nu);
            end

    end

else 

    error('Error: invalid simulation method \n');

end

% Store output
output.simM                                     = simM;
output.simX                                     = simX;
output.M_monthly                                = M_monthly;
output.X_monthly                                = X_monthly;
output.M                                        = Moms;

% Plot simulation output
if nargin>3 && to_plot==1
figure(fignum)
subplot(211)
plot(tvec,reshape(mean(simM,2),[Nt+1,Nu.Nreps])); hold on; 
plot(tvec_eom,reshape(mean(M_monthly,2),[size(M_monthly,1),Nu.Nreps]),'x'); hold on; drawnow;
title('M');
subplot(212)
plot(tvec,reshape(mean(simX,2),[Nt+1,Nu.Nreps])); hold on; 
plot(tvec_eom,reshape(mean(X_monthly,2),[size(M_monthly,1),Nu.Nreps]),'x'); hold on; drawnow;   
title('X');
end

end

